Testing Controls

How to interface with TCGA https://bioconductor.org/packages/devel/bioc/vignettes/TCGAbiolinks/inst/doc/index.html https://bioconductor.org/packages/release/workflows/vignettes/TCGAWorkflow/inst/doc/TCGAWorkflow.html#Environment

Loading the annotation file:

data("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
load("/Users/adiallo/Desktop/Dartmouth/Christiansen_lab/CF_Project/Methylation/Annotation/EPIC.hg19.manifest.RDATA")
attach(annotation)
manifest <- getManifest(IlluminaHumanMethylationEPICmanifest)

Loading the data and metadata information

targets <- read.csv("/Users/adiallo/Desktop/Thesis/Data_Documents/dm_57_samples.csv")
targets$patient <- paste(targets$Sentrix_ID,targets$Sentrix_Position,sep="_")
rownames(targets) <- targets$patient
targets$SampleID<- targets$patient

targets_expanded <- read.csv("/Users/adiallo/Desktop/Thesis/Data_Documents/data_all.csv")
targets_expanded$patient <- paste(targets$Sentrix_ID,targets$Sentrix_Position,sep="_")
rownames(targets_expanded) <- targets_expanded$patient
targets_expanded$SampleID<- targets_expanded$patient
targets = targets_expanded

idat = "/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/no_match/"
RGset_32 = read.metharray.exp("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/no_match/idats/",recursive = TRUE) 
#RGset_25 = openSesame(idat , func = getBetas) 
RGset_25 = read.metharray.exp("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/no_match/25_samples/",recursive = TRUE) 

Normalize my Data

Noob_25_m = preprocessNoob(RGset_25)
Noob_32_m = preprocessNoob(RGset_32)

Getting the beta values

Betas_25<-getBeta(Noob_25_m)
Betas_32<-getBeta(Noob_32_m)

Betas_25<- sesame::betasCollapseToPfx(Betas_25)
#Betas_32<- sesame::betasCollapseToPfx(Betas_32)
colnames(Betas_25) = colnames(Noob_25_m)

DH_CRC_Betas <- merge(Betas_25, Betas_32, by = "row.names")
rownames(DH_CRC_Betas) <- DH_CRC_Betas$Row.names
DH_CRC_Betas <- DH_CRC_Betas[,-1] 

Running EpiDish to obtain Fibroblast information

# Define the beta matrix and reference matrix
beta_matrix <- DH_CRC_Betas  # Your dataset
ref_matrix <- centEpiFibIC.m  # Reference matrix for EpiDISH (e.g., centEpiFibIC.m)

# Run EpiDISH
epidish_results <- epidish(beta.m = beta_matrix, ref.m = ref_matrix, method = "RPC")

# Extract the estimated cell type proportions
cell_type_proportions <- epidish_results$estF

# Convert the matrix to a data frame and scale proportions to percentages
cell_type_proportions_df <- as.data.frame(cell_type_proportions) * 100

# Add SampleID as a column from the row names
cell_type_proportions_df$SampleID <- rownames(cell_type_proportions)

# Inspect the results
head(cell_type_proportions_df)

Running HiTIMED to generate cell type proportions

HiTIMED_result<-HiTIMED_deconvolution(DH_CRC_Betas,"COAD",5,"tumor")

snapshotDate(): 2024-04-29
see ?HiTIMED and browseVignettes('HiTIMED') for documentation
loading from cache

HiTIMED_result$SampleID <- rownames(HiTIMED_result)

HiTIMED_result_immune<-HiTIMED_deconvolution(DH_CRC_Betas,"COAD",2,"tumor")

snapshotDate(): 2024-04-29
see ?HiTIMED and browseVignettes('HiTIMED') for documentation
loading from cache

HiTIMED_result_immune$SampleID <- rownames(HiTIMED_result)

Get GTEx data

Betas_GTEx = readRDS("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/GTEx/GTEx_samples.rds")
targets_GTEx <- read.csv("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/GTEx/targets.csv")
targets_GTEx$patient <- targets_GTEx$Filename
rownames(targets_GTEx) <- targets_GTEx$patient
targets_GTEx$SampleID<- targets_GTEx$patient

# --- Harmonize CpGs ---
common_cpgs_GTEx <- intersect(rownames(Betas_GTEx), rownames(centEpiFibIC.m))
betas_GTEx_filtered <- Betas_GTEx[common_cpgs_GTEx, ]
ref_filtered_GTEx <- centEpiFibIC.m[common_cpgs_GTEx, ]

# --- Run EpiDISH ---
epidish_GTEx <- epidish(beta.m = betas_GTEx_filtered, ref.m = ref_filtered_GTEx, method = "RPC")
cell_props_GTEx <- as.data.frame(epidish_GTEx$estF) * 100
cell_props_GTEx$SampleID <- rownames(cell_props_GTEx)

# Optional: Merge with metadata
cell_props_GTEx <- merge(cell_props_GTEx, targets_GTEx, by = "SampleID")

HiTIMED_result_GTEx<-HiTIMED_deconvolution(Betas_GTEx,"COAD",h = 5,"tumor")

snapshotDate(): 2024-04-29
see ?HiTIMED and browseVignettes('HiTIMED') for documentation
loading from cache

HiTIMED_result_GTEx$SampleID <- rownames(HiTIMED_result_GTEx)

HiTIMED_result_immune_GTEx<-HiTIMED_deconvolution(Betas_GTEx,"COAD",h = 2,"tumor")

snapshotDate(): 2024-04-29
see ?HiTIMED and browseVignettes('HiTIMED') for documentation
loading from cache

HiTIMED_result_immune_GTEx$SampleID <- rownames(HiTIMED_result_GTEx)

Get TCGA data

tcga.data <- readRDS("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/TCGA/TCGA_data.rds")
tcga.dnam <- tcga.data[["dnam"]]
tcga.pheno <- tcga.data[["pheno"]]
tcga.cell <- tcga.data[["cell_types"]]
msi_dat<-readRDS("/Users/adiallo/Desktop/Thesis/Data_Documents/All_Data/DNA_Methylation/dm_data/TCGA/msi_tcga.rds")
tcga.cell$COAD_Deconv1#$MSS

tcga.cell$COAD_Deconv2#MASS

tcga_T = t(tcga.dnam)

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggpubr)

# ---- Function to Reshape Data ----
extract_subset <- function(df, source_label, cell_types = c("CD8T", "Tumor", "Epithelial")) {
  df %>%
    select(all_of(cell_types), SampleID) %>%
    pivot_longer(-SampleID, names_to = "CellType", values_to = "Proportion") %>%
    mutate(Source = source_label)
}

# ---- Combine Datasets ----
df_DH <- extract_subset(HiTIMED_result, "DH")
df_GTEx <- extract_subset(HiTIMED_result_GTEx, "GTEx")
df_TCGA <- extract_subset(HiTIMED_result_TCGA, "TCGA")

combined_df <- bind_rows(df_DH, df_GTEx, df_TCGA)

# ---- Define comparisons ----
comparisons <- list(c("DH", "GTEx"), c("DH", "TCGA"))

# ---- Boxplot with Significance Annotations ----
ggplot(combined_df, aes(x = Source, y = Proportion, fill = Source)) +
  geom_boxplot(outlier.shape = NA, width = 0.7) +
  facet_wrap(~ CellType, scales = "free_y") +
  stat_compare_means(comparisons = comparisons, method = "wilcox.test",
                     label = "p.signif", hide.ns = TRUE) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Cell Type Proportions Across Datasets",
       y = "Estimated Proportion",
       x = "Dataset") +
  theme_minimal(base_size = 14) +
  theme(strip.text = element_text(size = 14),
        axis.text.x = element_text(angle = 30, hjust = 1))

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggpubr)

# ---- Function to Extract Cell Types ----
extract_subset <- function(df, source_label, cell_types = c("CD8T", "Tumor", "Epithelial")) {
  df %>%
    select(all_of(cell_types), SampleID) %>%
    pivot_longer(-SampleID, names_to = "CellType", values_to = "Proportion") %>%
    mutate(Source = source_label)
}

# ---- Combine Only DH and GTEx ----
df_DH <- extract_subset(HiTIMED_result, "DH")
df_GTEx <- extract_subset(HiTIMED_result_GTEx, "GTEx")

combined_df <- bind_rows(df_DH, df_GTEx)

# ---- Define Comparisons ----
comparisons <- list(c("DH", "GTEx"))

# ---- Plot ----
ggplot(combined_df, aes(x = Source, y = Proportion, fill = Source)) +
  geom_boxplot(outlier.shape = NA, width = 0.7) +
  facet_wrap(~ CellType, scales = "free_y") +
  stat_compare_means(comparisons = comparisons, method = "wilcox.test",
                     label = "p.signif", hide.ns = TRUE) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Comparison of Cell Type Proportions: DH vs GTEx",
       y = "Estimated Proportion",
       x = "Dataset") +
  theme_minimal(base_size = 14) +
  theme(strip.text = element_text(size = 14),
        axis.text.x = element_text(angle = 30, hjust = 1))

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggpubr)

# ---- Function to Extract HiTIMED Cell Types ----
extract_subset_hitimed <- function(df, source_label, cell_types = c("CD8T", "Tumor", "Epithelial")) {
  df %>%
    select(all_of(cell_types), SampleID) %>%
    pivot_longer(-SampleID, names_to = "CellType", values_to = "Proportion") %>%
    mutate(Source = source_label)
}

# ---- Extract HiTIMED Results for DH and GTEx ----
df_DH_hitimed <- extract_subset_hitimed(HiTIMED_result, "DH")
df_GTEx_hitimed <- extract_subset_hitimed(HiTIMED_result_GTEx, "GTEx")

# ---- Function to Extract Fibroblast Estimates from EpiDISH/EpiScore ----
extract_subset_fib <- function(epi_df, source_label) {
  epi_df %>% 
    select(Fib, SampleID) %>% 
    pivot_longer(-SampleID, names_to = "CellType", values_to = "Proportion") %>% 
    mutate(Source = source_label)
}

# ---- Extract Fibroblast Results for DH and GTEx ----
df_DH_fib <- extract_subset_fib(cell_type_proportions_df, "DH")
df_GTEx_fib <- extract_subset_fib(cell_props_GTEx, "GTEx")

# ---- Combine All Results ----
combined_df <- bind_rows(df_DH_hitimed, df_GTEx_hitimed, df_DH_fib, df_GTEx_fib)

# ---- Define Comparisons (DH vs GTEx for each cell type) ----
comparisons <- list(c("DH", "GTEx"))

# ---- Plot: Faceted Boxplots by CellType ----
ggplot(combined_df, aes(x = Source, y = Proportion, fill = Source)) +
  geom_boxplot(outlier.shape = NA, width = 0.7) +
  facet_wrap(~ CellType, scales = "free_y") +
  stat_compare_means(comparisons = comparisons, method = "wilcox.test",
                     label = "p.signif", hide.ns = TRUE) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Comparison of Cell Type Proportions: DH vs GTEx",
       y = "Estimated Proportion (%)",
       x = "Dataset") +
  theme_minimal(base_size = 14) +
  theme(strip.text = element_text(size = 14),
        axis.text.x = element_text(angle = 30, hjust = 1))

For obtaining more data from TCGA

query_meth <- GDCquery(
    project = "TCGA-COAD",
    data.category = "DNA Methylation",
    platform = "Illumina Human Methylation 450"
)

--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg38
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-COAD
--------------------
oo Filtering results
--------------------
ooo By platform
----------------
oo Checking data
----------------
ooo Checking if there are duplicated cases
Warning: There are more than one file for the same case. Please verify query results. You can use the command View(getResults(query)) in rstudio
ooo Checking if there are results for the query
-------------------
o Preparing output
-------------------

GDCdownload(query_meth)

data_type in query
Masked Intensities
Methylation Beta Value

Error in GDCdownload(query_meth) : 
  We can only download one data type. Please use data.type argument in GDCquery to filter results.